Past and future gauge in numerical relativity 
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ABSTRACT 



Numerical relativity describes a discrete initial value problem for general relativ- 
D ■ ity. A choice of gauge involves slicing space-time into space-like hyper surf aces. This 

introduces past and future gauge relative to the hypersurface of present time. Here, we 
propose solving the discretized Einstein equations with a choice of gauge in the future 
and a dynamical gauge in the past. The method is illustrated on a polarized Gowdy 
wave. 
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The initial value problem for general relativity is receiving much attention in the prediction of 
wave- forms from candidate sources for the upcoming gravitational wave detectors LIGO/VIRGO 
(1; 2). The structure of gravitational waves has recently been elucidated in new hyperbolic formula- 
tions of general relativity (see, e.g., (3) for references), which holds promise for accurate integration 
schemes. Long-time integrations also require accurate conservation of the associated elliptic con- 
straints, representing energy and momentum conservation (4). While analytically these constraints 
are conserved under dynamical evolution, the nonlinear nature of the equation typically tends to 
introduce some numerical departure from exact conservation. 

Gauge choice in numerical relativity involves slicing space-time into space-like hyper surf aces. 
Data from one hypersurface are evolved numerically onto the next, for example using a hyperbolic 
formulation. Algebraic slicing is particularly illustrative, which considers prescribed lapse and shift 
functions N a = gt a for the metric tensor g a \,. 

Here, we propose an new approach for numerical relativity with preservation of the constraints 
by through choice of gauge in the future and a dynamical gauge in the past. With prescribed gauge 
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and dynamical space-like components of the metric in the future, this obtains a complete system 
of ten equations in ten unknowns. 

The presented approach recognizes two types of constraints in Hamiltonian approaches within 
the context of the underlying four-covariant theory of general relativity. Recall that Hamilto- 
nian formulations assign the three- metric hij(t) and its canonical momentum iTij(t) as dynamical 
variables to a hypersurface of constant coordinate time t. This implicitly carries a geometric 
constraint, in having hij(t) as exact projections of the four-covariant metric g a b onto the hyper- 
surfaces of constant t. In addition, general relativity obeys energy-momentum conservations by 
the Bianchi identity, which translates into the Gauss-Codacci relations on the same hyper surfaces. 
These observations may be contrasted with the nonlinear wave equations of (6), which describe the 
four-covariant hyperbolic evolution of the Riemann-Cartan connections uj^ first, with subsequent 
slicing of space-time for the purpose of numerical implementation. In the continuum limit, both the 
projections and energy- momentum constraints are exact and consistent. This need not carry over 
to the discrete case, given the nonlinear nature of the Einstein equations. This may be illustrated 
by the observed violations of energy- momentum constraints in present numerical simulations (4) . 

A discretization of the Einstein equations introduces finite deviations from the continuum limit. 
For the purpose of numerical relativity, these deviations are acceptable, provided they permit stable 
evolution and convergence upon refinement of the discretization. The question therefore is: where 
are discretizations allowed to introduce deviations from the continuum theory? The above suggests 
to consider seeking a trade off between the geometric constraint of exact projections and energy- 
momentum conservation. This points towards a numerical scheme in which the past slicing is 
adjusted within the chosen discretization level, so as to satisfy energy- momentum conservation 
identically for a given definition of the discretized Einstein tensor. This obtains a complete system 
of evolution equations, as follows from simple counting: the contracted Bianchi identity shows that 
the Einstein equations provide six constraints on the 2 nd time-derivatives of the 3-metric hij and 
four constraints on the first time-derivatives of the lapse and shift functions N a = g a f. The first 
tells us that the discretized Einstein equations live on three time-slices, whereas the latter indicates 
that the lapse and shift functions are constraint on a pair of them, e.g., the first (past) and the 
third (future) time-slice. Thus, choosing a gauge for one of these leaves the gauge on the other 
determined self-consistently by the Einstein equations. It seems natural to maintain control over 
the future gauge, i.e., the future gauge is by choice of the the user. This leaves the past gauge as a 
dynamical variable, to be determined implicitly during numerical evolution. This can be achieved 
using Newton's method. 
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Non-exact projections naturally permit an uncertainty between the three-metric and its canon- 
ical momentum within the underlying context of a four-covariant theory, i.e.: also in regards to 
the association with the hypersurface at hand. In the covariant approach of (6), this would thus 
reflect an uncertainty in the tetrad elements, which define the projection, and their connections. 
This points towards a potential connection to quantum gravity. Indeed, soon after this work was 
proposed (5), the author learned of a very interesting independent discussion on the problem of 
consistent discretizations in this context (7). 



2. A discretized initial value problem 

We illustrate our this approach on the vacuum Einstein equations, described by the vanishing 
Ricci tensor 

Rab = 0. (1) 

The Ricci tensor R a b is a second-order expression in the metric g a b- Hence, (1) defines a relationship 
between metric data {g™ h _1 , g™ b , fl^ 1 ) on a triple of time-slices i n _i < t n < t n+ i: 

R ab (9 n ab +1 ,g:b,9 n ab T 1 )=0. (2) 

Here, Rbd = R a bed derived from the Riemann tensor 

R a bed = ddFbc - dc^td + r^Ted - rf d rg C . (3) 

This expression (3) can be discretized by finite differencing on a triple of time-slices with preserva- 
tion of the quasi-linear second-order structure of R a b- 

Algebraic gauge-fixing takes the form of specifying the components N a = gt a in coordinates 
{x a }a = i with t = x 1 time-like. A gauge-choice on a triple of time-slices, therefore, amounts to 
a choice of (N^ -1 , N™, N™ +1 ). Recall that this gauge-choice in the metric arises explicitly in the 
Gauss-Codacci relations for energy-momentum conservation. The components hij = gij, where 
i,j = 2,3,4 refer to projections of the metric into the time-slice t =const., which describe the 
dynamical part of the metric. The combination (hij,N a ) reflects the mixed hyperbolic-elliptic 
structure in numerical relativity and (1) represents ten evolution equations in these variables on a 
triple of time-slices. 

In algebraic gauge-fixing, we prescribe N™ +1 as a future gauge in computing h^ 1 on a future 
hypersurface t = t n+ \ from data at present and past hypersurfaces t = t n -\ and t = t n . We propose 
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closure by re-introducing iV™" 1 as dynamical gauge in the past, leaving fixed. Combined, this 
defines an advanced hyperbolic-retarded elliptic evolution of the metric. The paritioning of the 
metric in past and future variables as 
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thus obtains ten dynamical variables in the ten equations 

Rab(K+ 1 ,K-\---) = at t = t n . (5) 

Here the dots refer to the remaining data (h^ 1 , hfj , TV™ , iV" +1 ), which are kept fixed while solving 
for {h™* 1 , iV" _1 ). Thus, (5) which takes into account all ten Einstein equations with no reduction of 
variables. Time-stepping by (5) evolves the metric into the future with dynamical gauge in the past, 
in an effort to satisfy energy-momentum conservation within the definition of the discretized Ricci 
tensor. Because (5) comprises derivatives of N a only to first-order in time, numerically though the 
data -/V™ +1 and A r ™ _1 , we anticipate that the evolution of N a is of first-order in the t— discretization 
At. This introduces non-exactness in /i™ -1 as projections of g a f, on t — At to within the same order 
of accuracy. It may result in a first-order drift in the t— labeling of the hypersurfaces - permitted 
by coordinate invariance. 

The presented approach can be illustrated on a polarized Gowdy wave. Gowdy cosmologies 
are an extensively studied class of universes with compact space-like hypersurfaces with two Killing 
vectors d a and d$. With cyclic boundary conditions, the space-like hypersurfaces are homeomorphic 
to the three-torus as a model universe collapsing into a singularity. The associated line-element is 
(see, e.g., (8)) 

ds 2 = e (r-A)/2 (_ e ~2r dr 2 + ^2) + ^ (6) 

where A = A(r, 6) and gE denotes the surface element in the space supported by the two Killing 
vectors. Polarized Gowdy waves form a special case, which permit a reduction to 

dT? = e- T (eW + e' p d5 2 ) . (7) 

Here P satisfies a linear wave-equation P TT = e~ 2r Pqq\ a long wave-length solution is 



PoM) = Y o (e- T )cos0, 



(8) 
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where Yq is the Bessel function of the second kind of order zero. This leaves 

A(r, 6) = ±Yo(e^)Yi(e- r )e- T cos 29 + ± jf ^ (Y^(s) + Y \s)) sds. (9) 

A spectrally accurate numerical integration is described in (9). 

The implicit equation (5) for the dynamical variables (h™* 1 , N£ _1 ) has been implemented 
numerically. We have done so by solving for the all ten components (h™^ 1 , iV™ -1 ) using Newton 
iterations on these variables. This procudere uses a numerical evaluation of the Jacobian 

- H « io > 

where the capital indices A, B = 1, 2, • • • , 10 refer to the labeling 

Ra = {Rll, R22, #33) R44, Rl2, Rl3, Rl4, R23, R24, R34), /^x 

U B = (iVr 1 , <\ ^ 2 +1 , h^ 1 ,!^- 1 ,^- 1 , N2~\ h n 2 $\ h£\ Kit 1 ). 1 j 

The Ricci tensor (3) has been implemented by second-order finite differencing, such that it remains 
quasi-linear in the second derivatives. In particular, the Christoffel symbols 

r afe = -j9 Ce {9cb,a + 9ac,b ~ 9ab,c) (12) 

is obtained by symmetric finite-differencing on the metric components, and itself differentiated by 
the product rule following individual numerical differentiations of g ab and (g c b, a + 9ac,b ~ 9ab ,c)- The 
choice of future gauge iV™ +1 is provided by the the components 

9at = (e {T - X)/ \ 0,0.0) (13) 

of the analytical line-element (6-9), which facilitates error analysis by direct comparison of the 
numerical results with the analytic expression for the line-element (6). It will be appreciated that 
in principle other choices of iV™ +1 can be made. 

Fig. 1 shows numerical results for evolution of initial data on < r < 4. The results show that 
all Einstein equations in the form of R a b = are satisfied with arbitrary precision, while the metric 
components are solved accurately to within one percent. The asymptotic behavior of the implicit 
corrections to the lapse functions are shown in Fig. 2. Note that these corrections are finite to 
first-order in At, a testimonial to the appearance of the lapse function in the Einstein equations to 
first-order in time. 



In summary, a dynamical gauge in the past gives a complete number often degrees of freedom in 
evolving to a future hypersurface, permitting all ten Einstein equations to be satisfied with arbitrary 



-6- 



precision. The Einstein equations are hereby satisfied numerically on a triple of hypersurfaces of 
past, present and future within the definition of the discretized Riemann tensor. The simulation of 
a nonlinear one-dimensional Gowdy wave by implicit time-stepping according to the ten discretized 
vacuum Einstein equations (5) serves to illustrate a numerical implementation. The presented 
approach ensures conservation of energy and momentum, within the definition of the discretization 
used. Satisfying these constraints is generally a necessary condition for stability in long-time 
integrations. It is an important open question if satisfying these constraints is a sufficient condition 
for stability. It would be of interest to study this method in evolving Schwarzschild black holes 
with singularity-avoiding slicings of spacetime, and its behavior in the presence of outgoing radiative 
boundary conditions. More generally, it would be of interest to consider this approach in higher 
dimensions, including a self-consistent integration of any of the modern hyperbolic formulations 
and efficient elliptic solvers. 
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Figure Captions 

FIGURE 1. Shown is the simulation for < r < 4 of the polarized Gowdy wave. The solutions 
P(t,6) and X(t,9) are displayed as a function of (r, 6>) (upper windows). The middle windows 
display the solutions for r = 4, wherein the circles denote the numerical solution and the solid lines 
the analytical solution. The r— evolution of the errors (lower windows) are computed relative to 
the analytical solution to Gowdy's reduced wave equation. The simulations discretize 9 by m\ = 64 
points and the r— interval by mi = 1024 time-steps. Particular to the proposed numerical algorithm 
is a dynamical gauge in the past and a prescribed gauge in the future. This permits satisfying all of 
the discretized Einstein equations R a b = to within arbitrary precision by Newton iterations. The 
slight increase in the error of about 10 -10 reflects the exponential growth of the analytic solution, 
because the Gowdy cosmology evolves towards a singularity. 

FIGURE 2. Shown are the self-consistent corrections on the slicing t = t n+ ±, introduced by the 
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difference between the past gauge iV" - " 1 (i n +2) to the hypersurface t = t„ + 2 and the earlier future 
gauge iV n+1 (i n ) to the hypersurface t = t n . The three curves refer to different discretizations mi = 
16, 32 and 64 points with, respectively, 1712 = 256, 512 and 1024 time-steps. These similar results for 
various discretizations indicate asymptotic behavior consistent with the first-order appearance of the 
lapse function in the Einstein equations. A first-order accuracy in lapse introduces a commensurate 
offset in slicing or, equivalently, an offset in the coordinate t. 
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FIGURE 1 
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J n+1 (x n )-N n - 1 (x n+2 )|/N n (x n+1 )Ax 




FIGURE 2 



